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1  Summary  of  accomplishments 

The  overarching  goal  of  the  originally  proposed  work  was  the  development  of  a  new 
generation  scalable  Lagrangian  finite  element  software  fi-amework  capable  of  capturing 
local  and  non-local  fracture  processes  in  solids  and  structures.  The  principal  objectives 
of  the  work,  as  stated  in  the  proposal,  were: 

L  Development  a  scalable  Lagrangian  software  fiamework  to  simulate  fiacture  and 
fragmentation  due  to  weapon-target  interactions. 

2.  Implementation  of  complex  phenomenological  as  well  as  physics-based  advanced 
material  strength  and  damage  models, 

3.  Validation  and  verification  of  software  through  simulation  of  several  benchmark 
DoD  application  problems 

All  the  objectives  have  been  successfully  achieved.  Some  of  the  main  developments 
are  summarized  in  the  text. 


2  The  scalable  Lagrangian  finite  element  software  frame¬ 
work 

Over  the  past  decade,  the  PI  has  developed  state-of-the-art  Lagrangian  finite  element 
(FE)  software  capable  of  simulating  crack  initiation  and  growth  in  solids  and  structures. 
The  software,  written  In  high-level  languages  (Fortran,  C  and  C++),  was  designed  to 
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execute  reliably  in  imi-processor  computing  environments  (personal  computers,  w'ork- 
stations,  etc.).  In  addition,  over  the  years,  it  underwent  thorou^  testing  and  perfor¬ 
mance  tuning.  Therefore,  diis  serial  software  become  a  suitable  candidate  to  be  turned 
into  a  core  of  the  parallel  FE  fiamewoik. 

Our  parallel  implementation  was  based  on  a  dynamic  decomposition  of  the  compu¬ 
tational  domain  among  the  processor  with  standard  partitioning  techniques  (METIS). 
A  two  level  approach  was  adopted  in  order  to  achieve  scalability  in  the  distributed 
computing  environments.  Initially,  a  set  of  multi-threaded  (MT)  extensions  to  the  se¬ 
rial  software  was  developed  as  to  provide  for  the  optimal  performance  on  symmetric 
multi-processors  (SMP).  Subsequently,  the  implementation  was  amended  by  the  ad¬ 
dition  of  extensions  based  on  the  standard  message  passing  interface  (MPI),  designed 
to  allow  for  the  distribution  of  simultaneous  computational  tasks  across  heterogenous 
parallel  computing  environments.  This  two  tier  approach  lead  to  exceptional  scalability 
and  performance  (Figure  1). 


Figure  1:  Measured  speedup  of  the  pardlel  software  fmmework  for  the  varying  prob¬ 
lem  size. 

In  the  design  process,  we  decided  to  consider  the  use  of  commercial  products  for 
the  pre-  and  post-processing.  To  this  end  we  developed  an  interface  for  the  FEMAP 
neutral  file  format  as  the  model  input.  FEMAP  is  a  popular  engineering  FE  analysis 
package  that  supports  solid  modeling,  FE  model  generation  and  post-processing.  In 
addition,  FEMAP  offers  the  ability  to  import  pre-existing  FE  models  in  many  CAD  file 
formats.  In  an  attempt  to  allow  the  end-user  to  create  high-quality  visualizations,  an 
interface  for  the  output  in  TECPLOT  visualization  toolkit  file  format  was  also  built. 

The  material  library  is  an  integral  and  essential  part  of  the  parallel  FE  fi*amework. 
In  order  to  facilitate  the  reuse  or  exchange  existing  material  models,  in  particular  those 
models  that  had  been  implemented  in  other  DoD  software,  the  strict  material  inter¬ 
face  guidelines  (MEG)  were  established  and  followed.  This  MIG  format  allowed  the 
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US  Army  Research  Laboratory  (ARL)  personnel  to  successfully  import  the  material 
models  implemented  in  other  FE  software  (e.g.  EPIC), 

The  software  development  required  close  collaboration  between  both  the  Caltech 
and  the  ARL  teams.  In  order  to  simplify  the  process  a  CVS  source  repository  was 
set  up,  with  the  ARL  team  having  full  access  to  the  parallel  framework  source  code 
through  out  the  entire  course  of  the  project. 

3  Parallel  fracture  algorithm 

Many  engineering  problems  involving  fracture  and  fragmentation  phenomena  cannot 
be  addressed  on  serial  and  vector  computers  due  to  processing  time  or  storage  require¬ 
ments.  However,  in  the  context  of  MPI-based  domain  decomposition  algorithm,  a  sim¬ 
ple  extension  of  the  serial  fracture  and  fragmentation  approach  is  not  feasible.  There¬ 
fore,  as  an  initial  step  towards  the  scalable  fully  three-dimensional  (3D)  fracture  and 
fragmentation  software  implementation,  a  considerably  simpler,  two-dimensional  (2D) 
version  was  developed.  On  the  basis  of  this  effort,  we  adopted  an  approach  based  on 
the  introduction  of  cohesive  elements  in  the  FE  discretization  directly  in  the  vicinity 
of  all  processor  boundaries  at  the  onset  of  the  simulation.  This  approach  significantly 
reduced  the  complexity  of  the  3D  fracture  and  fragmentation  software  module,  as  well 
as  decreased  the  development  time.  In  addition,  the  3D  implementation  was  designed 
as  easily  expandable  as  to  allow  to  incorporate  other  models  of  physical  behavior  that 
requires  topological  changes  in  the  FE  discretizations  (e.g,  adiabatic  shear  bands). 


4  Adiabatic  shear-band  modeling 

Solids  deforming  at  high  rates  often  develop  narrow  layers  of  intense  shearing.  Out¬ 
standing  features  of  these  dynamic  shear  bands  are  their  thinness,  high  local  shear 
strains,  ultra-high  local  shear  strain  rates,  significant  local  temperature  rises  and  high 
propagation  speeds.  When  simulating  shear-bands,  the  commonly  adopted  approach 
calls  for  mesh  adaption  to  resolve  the  fine  scales  in  the  solution.  This,  in  turn,  may 
significantly  increase  the  size  of  the  FE  problem,  placing  more  demands  in  terms  of 
computer  power  or  storage.  An  alternative  approach  is  to  employ  a  special  type  of  finite 
elements,  which  serve  the  purpose  of  capturing  sub-grid  localization  processes  such  as 
shear  bands  and  void  sheets[8].  These  elements  consist  of  two  surfaces,  attached  to 
the  abutting  volume  elements,  which  can  separate  and  slip  relative  to  each  other.  The 
kinematics  of  the  strain-localization  elements  is  three-dimensional  and  allows  for  finite 
deformations.  The  behavior  of  the  strain-localization  elements  is  governed  directly  by 
the  same  constitutive  relation  which  describes  the  deformation  of  the  volume  elements. 
In  particular  the  localization  elements  allow  for  the  development  of  displacement  and 
velocity  jumps  across  volume  element  boundaries.  The  thickness  of  the  localized  zone 
is  set  by  an  additional  field  variable  which  is  determined  variationally.  The  localization 
elements  are  inserted,  and  become  active,  only  when  localized  deformations  become 
energetically  favorable. 
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5  Tetrahedral  composite  finite  elements 

The  quadratic  tetrahedral  elements  have  been  know  to  under-perform  in  explicit  time 
integration  and  contact-impact  problems,  due  to  the  comer  node  masses  being  negative 
or  zero,  Thoutireddy  et  aL  [6]  introduced  a  composite  ’CT3D’  tetrahedral  element  that 
overcomes  the  difficulties  inherent  to  the  quadratic  elements,  making  it  perfectly  suited 
for  impact  problems.  The  composite  tetrahedron  consists  of  an  ensemble  of  12  four- 
node  linear  tetrahedral  elements,  coupled  to  a  linear  assumed  deformation  defined  over 
the  entire  domain  of  the  composite  element.  The  formulation  of  the  element  accounts 
for  hilly  non-linear  kinematics.  Integrals  over  the  domain  of  the  element  are  computed 
by  a  five-point  quadrature  rule.  The  element  passes  the  patch  test  in  arbitrarily  distorted 
configurations.  Numerical  tests  demonstrate  that  CT  element  has  been  found  to  possess 
a  convergence  rate  comparable  to  those  of  linear  simplicial  elements,  and  that  these 
convergence  rates  are  maintained  as  the  near-incompressible  limit  is  approached. 


6  Validation  tests 

We  have  validated  the  new  three-dimensional  parallel  c^ability  by  means  of  extensive 
comparison  with  archival  data.  Some  relevant  work  in  this  direction  is  the  following: 

6,1  Three-dimensional  modeling  of  dynamic  mixed-mode  fracture 

As  the  first  validation  case,  we  have  taken  the  dynamic  mixed-mode  fracture  as  a  con¬ 
venient  basis  for  assessing  the  predictive  ability  of  cohesive  models  in  applications 
that  involve  competition  between  crack  initiation  and  growth  [5].  The  predictions  of 
our  numerical  model  are  validated  against  the  experimental  data  of  John  and  Shah  [4]. 
These  tests  were  concerned  with  three-point-bend  beam  (TPB)  plain  concrete  speci¬ 
mens  containing  a  pre-crack  shifted  firom  the  central  cross-section.  The  location  of 
the  pre-crack  was  characterized  by  the  relative  offset  7  —  2offset/L,  where  L  was 
the  length  of  the  beam.  The  presence  of  the  pre-crack  led  to  asymmetric  loading  con¬ 
ditions  and  resulted  in  mixed-mode  fracture.  Specimens  were  tested  dynamically  by 
means  of  a  modified  Chaipy  tester.  By  increasing  its  offset,  the  pre-crack  was  subjected 
to  increasingly  mode-II  fracture  conditions,  and  the  trajectory  of  the  crack  varied  ac¬ 
cordingly.  In  addition,  the  experimental  data  is  indicative  of  a  competition  between 
two  fracture  mechanisms:  the  growth  of  the  pre-crack;  and  the  nucleation  and  subse¬ 
quent  growth  of  a  crack  within  the  central  cross-section  of  the  specimen.  The  fracture 
behavior  of  the  specimen  was  also  observed  to  be  sensitively  dependent  on  the  rate  of 
loading. 

Figure  3  displays  computed  crack  patterns  for  relative  offsets  7  =  0.6, 0.7  and  0.8, 
respectively,  and  thus  illustrates  the  dependence  of  the  predicted  &ilure  mode  on  the 
location  of  the  notch.  In  the  experiments  of  John  as  Shah  [4]  a  change  in  &ilure  mode 
from  pre-notch  extension  to  crack  nucleation  within  the  plane  of  the  load  was  observed 
at  a  critical  offset  of  7  =  0.7.  Our  calculations  are  in  very  good  agreement  with  the 
experimental  observations,  as  the  computed  critical  7  is  identical  to  its  experimental 
counterpart.  Moreover,  our  results  show  that  the  model  is  indeed  capable  of  closely 
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Figure  2:  The  three-point-bend  beam  configuration  of  John  and  Shah  [4]  (dimensions 
in  mm). 


capturing  the  observed  fracture  patterns,  crack  trajectories  and  crack-tip  displacement 
fields  corresponding  to  different  experimental  configurations.  In  particular,  it  correctly 
accounts  for  the  relative  stability  regimes  of  the  competing  pre-crack  growth  vs.  crack 
nucleation  mechanisms. 

6.2  Ring  expansion  and  fragmentation 

The  second  validation  case  concerns  the  simulation  of  forced  expansion  and  ductile 
fracture  and  fragmentation  of  U-6%Nb  rings  [7,  1].  A  sample  ring  expands  driven 
by  a  driver  ring  under  the  magnetic  field  induced  by  the  electric  current  ^plied  to  a 
solenoid.  The  ring  has  an  inner  diameter  of 34.37  mm,  an  outer  diameter  of 35.89  mm, 
and  thickness  of  0,76  mm.  An  arrestor  ledge  in  the  experimental  device  prevents  the 
driver  ring  to  expand  beyond  a  diameter  of  40.0  mm.  In  contrast,  the  sample  ring  can 
expand  fi*eely,  until  it  eventually  fragments.  The  initial  velocity  history  ^plied  to  the 
U-6%Nb  ring  varies  according  to  the  three  different  accelerating  voltages  of  5.0  kV, 
6.0  kV,  and  7.5kV[l,3]„ 

The  U-6%Nb  alloy  is  assumed  to  obey  the  porous  plasticity  model  due  to  Weinberg 
et  al  [7],  with  the  strain  localization  part  of  Yang  et  al  [8].  Figure  4  shows  the  finite 
element  discretization  employed  in  our  calculations.  The  initial  number  of  nodes  is 
2, 634  and  the  number  of  elements  is  882.  The  nucleation  and  propagation  of  fi’acture 
is  simulated  by  adaptively  inserting  strain  localization  elements  when  the  accumulated 
plastic  strain  in  the  bulk  elements  reaches  the  value  of  0.4, 

The  experimental  velocity  histories  reported  by  Becker  [1],  as  well  as  those  ob¬ 
tained  in  our  simulations,  are  shown  in  Figure  5,  In  order  to  simulate  the  effect  of  the 
driver  ring,  the  radial  velocity  of  all  the  nodes  located  on  the  iimer  radius  is  initially  pre¬ 
scribed  to  match  the  radial  velocity  of  the  driver  ring,  as  measured  in  the  experiment, 
up  to  the  peak  velocity  for  each  curve  in  Figure  5.  A  noticeable  decrease  in  the  ring 
velocity  may  be  observed,  a  clear  symptom  of  the  kinetic  energy  being  converted  into 
the  work  of  plastic  deformation.  This  regime  is  subsequently  followed  by  a  plateau, 
indicative  of  ring  fi'agmentation.  Clearly,  our  calculations  are  capable  of  capturing 
the  onset  of  the  plastic  flow,  as  well  as,  the  initiation  of  fracture  and  fragmentation. 
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Figure  3:  Numerical  simulations  of  the  three-point-bend  experiments  of  John  and 
Shah  [4],  EfiFect  of  the  pre-notch  offset  on  crack  pattern  and  feilure  mode.  The  color 
coding  corresponds  to  the  magnitude  of  the  vertici  displacement. 
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Figure  4:  Finite  element  mesh  for  U-6%Nb  rings:  2, 634  nodes,  882  elements. 


-inn 


Time  (s) 


Figure  5:  Velocity  history  for  U-6%Nb  rings.  Experimental  data  taken  from  Becker 
[1]. 
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Moreover,  our  results  are  in  excellent  agreement  with  the  finite  element  simulations  of 
Becker  [1],  despite  the  considerably  coarser  discretization  of  the  ring  (882  vs.  15, 000 
elements). 


Figure  6:  Dynamic  fracture  and  fragmentation  of  U-6%Nb  rings.  Configurations  at 
OfjLS,  30fj.s,  60iJ>s  and  90fis. 


A  sequence  of  frames  of  the  fracture  and  fragmentation  of  the  U-6%Nb  ring  for  the 
7,5kV  applied  voltage  is  shown  in  Figure  6.  The  time  elapsed  between  fiames  is  30 
and  the  color  bar  is  related  to  the  the  magnitude  of  the  accumulated  plastic  strain. 
For  this  particular  voltage,  the  time  to  fracture  is  approximately  30  /us  [1].  The  figure 
clearly  shows  the  e?q)ansion  of  the  ring  and  the  process  of  fiacture  and  fragmentation 
that  begins  around  30  ^s.  The  number  of  major  fragments  obtained  in  the  simulation 
is  25,  in  good  agreement  with  the  total  of  19  fragments  recovered  by  Becker  [1]. 

6.3  Ballistic  Penetration  of  Mild-Steel  Plate  -  Plugging 

The  final  validation  case  involves  the  formation  of  a  shear  plug  in  a  circular  steel  plate 
subjected  to  an  impact  by  a  blunt  -nosed  projectile.  In  the  experiments  of  Borvik  et 
al  [2],  the  projectile  launched  from  a  gas  gun  impacted  the  center  of  the  target  steel 
plate  at  296  m/s.  The  plate  was  clamped  in  a  rigid  circular  frame.  The  recorded  data 
included  the  initial  and  residual  velocities  at  the  free  end  of  the  projectile. 

The  formation  of  shear-bands  plays  an  essential  role  in  the  penetration  of  the  target 
plate  by  the  projectile,  as  it  establishes  the  clear  shear  failure  path.  Should  the  cre¬ 
ation  of  the  shear-bands  in  the  model  be  constrained,  the  target  plate  does  not  become 
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perforated  as  the  plastic  deformation  is  distributed  over  larger  area. 

The  target  plate  has  the  diameter  of  500  mm  with  the  thickness  of  10  mm.  The 
length  of  the  projectile  is  80  mm  and  the  diameter  20  mm.  We  model  both  the  target  and 
projectile  material  as  J2“Plasticity  mild-steel  with  thermal  updates.  The  finite  element 


Figure  7:  The  finite  element  mesh  for  ballistic  penetration:  11, 594  nodes,  6, 430  ele¬ 
ments. 


mesh  used  in  our  simulations  is  shown  in  Figure  7.  The  initial  number  of  nodes  and 
elements  in  the  mesh  is  11, 594  and  6, 430,  respectively. 

The  con:q)uted  residual  projectile  velocity  at  38.8  /is  after  the  impact  is  251.3  m/s, 
in  contrast  with  the  value  of  217  m/s  reported  by  Borvik  et  al  [2].  We  attribute  the 
difference  to  the  feet  that  in  the  simulations  the  velocity  is  recorded  at  a  single  node 
location.  However,  the  average  over  a  set  of  nodes  yields  the  value  of  220.3  m/s,  now 
in  very  good  agreement  with  the  experimental  data.  In  Figure  8  we  show  a  sequence  of 
frames  during  the  perforation  of  the  target  plate  at  0/is,  33//s,  and  99fMS. 
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